Models for Q2
Erick Calderon-Morales and Leland Werden
library(dplyr)
library(nlme)
library(purrr)
library(performance)
library(tibble)
library(here)
library(reactablefmtr)
library(reactable)
library(emmeans)
library(car)
library(ggplot2)
library(modelr)
library(knitr)Load functions and data
# Load all joined dataset
source("./scripts/code_join_data_full_dataset.R")## [1] "Files loaded: data_rgr_cleaned"
## [1] "Files loaded: data_for_models"# Step was done like this because I am working with a subset of the data
# source cleaned data
source("./scripts/code_clean_data_nodules.R")## [1] "Files loaded: data_for_models" "Files loaded: data_nodules_cleaned"# Load custom made functions
source("./R/functions_models.R")
source("./R/function_plots.R")
source("./R/function_validation_plots.R")
source("./R/function_for_inference_emmeans_and_percentage_diff.R")Q2: How does increased nutrient and/or water availability influence seedling water- and nutrient-use traits and the relationships with N-fixing bacteria?
Model fitting
Traits
# Take response variables names
response_vars_q2 <-
set_names(c("amax", "gs", "sla"))## Create empty list
models_q2 <- list()model_q2_amax <- lme(amax ~ nfixer*treatment +init_height,
random = ~1|spcode,
data = data_for_models)
model_q2_amax <- list(model_q2_amax)
names(model_q2_amax) <- "amax"model_q2_gs <- lme(gs ~ nfixer*treatment +init_height,
random = ~1|spcode,
data = data_for_models)
model_q2_gs <- list(model_q2_gs)
names(model_q2_gs) <- "gs"model_q2_sla <- lme(sla ~ nfixer*treatment +init_height,
random = ~1|spcode,
data = data_for_models)
model_q2_sla <- list(model_q2_sla)
names(model_q2_sla) <- "sla" ## PNUE
model_q2_pnue_log <- lme(log(pnue) ~ nfixer*treatment +init_height,
random = ~1|spcode,
data = data_for_models)
model_q2_pnue_log <- list(model_q2_pnue_log)
names(model_q2_pnue_log) <- "pnue_log"model_q2_pnue_nlme <- lme(pnue ~ nfixer*treatment + init_height,
random = ~1|spcode,
weights = varIdent(form = ~1|spcode*treatment),
data = data_for_models)
model_q2_pnue_nlme <- list(model_q2_pnue_nlme)
names(model_q2_pnue_nlme) <- "pnue_nlme"## Narea_g_m2 log model
model_q2_n_area_log <- lme(log(narea_g_m2) ~ nfixer*treatment + init_height,
random = ~1|spcode,
data = data_for_models)
model_q2_n_area_log <- list(model_q2_n_area_log)
names(model_q2_n_area_log) <- "n_area_log"model_q2_n_area_nlme <- lme(narea_g_m2 ~ nfixer*treatment + init_height,
random = ~1|spcode,
weights = varIdent(form = ~1|spcode*treatment),
data = data_for_models)
model_q2_n_area_nlme <- list(model_q2_n_area_nlme)
names(model_q2_n_area_nlme) <- "n_area_nlme"## WUE log model
model_q2_wue_log <- lme(log(wue) ~ nfixer*treatment + init_height,
random = ~1|spcode,
data = data_for_models)
model_q2_wue_log <- list(model_q2_wue_log)
names(model_q2_wue_log) <- "wue_log"model_q2_wue_nlme <- lme(wue ~ nfixer*treatment + init_height,
random = ~1|spcode,
weights = varIdent(form = ~1|spcode*treatment),
data = data_for_models)
model_q2_wue_nlme <- list(model_q2_wue_nlme)
names(model_q2_wue_nlme) <- "wue_nlme"Nodule colonization data
# Delete unused variables
data_nodules_cleaned <-
data_nodules_cleaned %>%
# add id to rownames for keep track of the rows
column_to_rownames("id") %>%
dplyr::select(spcode, treatment, everything())Nodule weight
nlme_nodule_weight <- lme(estimated_total_nodule_mass_per_plant ~ treatment + init_height,
random = ~1|spcode,
weights = varIdent(form = ~1|spcode*treatment),
data = data_nodules_cleaned)
model_q2_nodule_weight <- list(nlme_nodule_weight)
names(model_q2_nodule_weight) <- "nodule_weight"Nodule count
nlme_nodule_count <- lme(log(total_number_of_plant_nodules) ~ treatment + init_height,
random = ~1|spcode,
weights = varIdent(form = ~1|spcode),
data = data_nodules_cleaned)
model_q2_nodule_count <- list(nlme_nodule_count)
names(model_q2_nodule_count) <- "nodule_count"Plant biomass investment to nodules
# Get Nfixers biomass
data_fixers_plant_biomass <-
data_for_models %>%
# Filter only Nfixers
filter(nfixer == 'fixer') %>%
# Select only the columns that will be used
select("spcode", "treatment", "init_height", "total_biomass") %>%
# Get plant's id
rownames_to_column("id")# Join datasets
data_plant_inves_to_nodule <-
data_nodules_cleaned %>%
# Get plant's id
rownames_to_column("id") %>%
# join data
full_join(., data_fixers_plant_biomass, by = join_by(id, spcode,
treatment,
init_height)) %>%
# Remove NA's
drop_na() %>%
# Calculate investment to nodules
mutate(nodule_mass_ratio = (estimated_total_nodule_mass_per_plant / total_biomass)*100)nlme_nodule_mass_ratio <- lme(log(nodule_mass_ratio) ~ treatment + init_height,
random = ~1|spcode,
weights = varIdent(form = ~1|spcode),
data = data_plant_inves_to_nodule)
model_q2_nodule_mass_ratio <- list(nlme_nodule_mass_ratio)
names(model_q2_nodule_mass_ratio) <- "nodule_mass_ratio"Append all models into a single list
# Append models to model list
models_q2 <- append(model_q2_amax, models_q2)
models_q2 <- append(model_q2_gs, models_q2)
models_q2 <- append(model_q2_sla, models_q2)
models_q2 <- append(model_q2_n_area_log, models_q2)
models_q2 <- append(model_q2_n_area_nlme, models_q2)
models_q2 <- append(model_q2_pnue_log, models_q2)
models_q2 <- append(model_q2_pnue_nlme, models_q2)
models_q2 <- append(model_q2_wue_log, models_q2)
models_q2 <- append(model_q2_wue_nlme, models_q2)
models_q2 <- append(model_q2_nodule_count, models_q2)
models_q2 <- append(model_q2_nodule_weight, models_q2)
models_q2 <- append(model_q2_nodule_mass_ratio, models_q2)names(models_q2)## [1] "nodule_mass_ratio" "nodule_weight" "nodule_count"
## [4] "wue_nlme" "wue_log" "pnue_nlme"
## [7] "pnue_log" "n_area_nlme" "n_area_log"
## [10] "sla" "gs" "amax"Model Assumptions
Maximal photosynthesis
validation_plots(models_q2$amax, data = data_for_models, group = "spcode")
Stomatal Conductance
validation_plots(models_q2$gs, data = data_for_models, group = "spcode")
SLA
validation_plots(models_q2$sla, data = data_for_models, group = "spcode")
Water Use Efficiency
validation_plots(models_q2$wue_log, data = data_for_models, group = "spcode")
validation_plots(models_q2$wue_nlme, data = data_for_models, group = "spcode")
PNUE
validation_plots(models_q2$pnue_log, data = data_for_models, group = "spcode")
validation_plots(models_q2$pnue_nlme, data = data_for_models, group = "spcode")
Nitrogen concentration per unit of area
validation_plots(models_q2$n_area_log, data = data_for_models, group = "spcode")
validation_plots(models_q2$n_area_nlme, data = data_for_models, group = "spcode")
Nodule weight
validation_plots(models_q2$nodule_weight, data = data_nodules_cleaned, group = "treatment")
Nodule count
validation_plots(models_q2$nodule_count, data = data_nodules_cleaned, group = 'spcode')
Plant invesment to nodules
validation_plots(models_q2$nodule_count, data = data_nodules_cleaned, group = 'spcode')
Model inference
## r2 models
models_q2 %>%
map(., r2_nakagawa) %>%
unlist()## nodule_mass_ratio.R2_conditional.Conditional R2
## 0.77835
## nodule_mass_ratio.R2_marginal.Marginal R2
## 0.09318
## nodule_weight.R2_conditional.Conditional R2
## 0.41047
## nodule_weight.R2_marginal.Marginal R2
## 0.13936
## nodule_count.R2_conditional.Conditional R2
## 0.99211
## nodule_count.R2_marginal.Marginal R2
## 0.05918
## wue_nlme.R2_conditional.Conditional R2
## 0.90104
## wue_nlme.R2_marginal.Marginal R2
## 0.83953
## wue_log.R2_conditional.Conditional R2
## 0.71998
## wue_log.R2_marginal.Marginal R2
## 0.55827
## pnue_nlme.R2_conditional.Conditional R2
## 0.54153
## pnue_nlme.R2_marginal.Marginal R2
## 0.25717
## pnue_log.R2_conditional.Conditional R2
## 0.49405
## pnue_log.R2_marginal.Marginal R2
## 0.19552
## n_area_nlme.R2_conditional.Conditional R2
## 0.85290
## n_area_nlme.R2_marginal.Marginal R2
## 0.50938
## n_area_log.R2_conditional.Conditional R2
## 0.86656
## n_area_log.R2_marginal.Marginal R2
## 0.58168
## sla.R2_conditional.Conditional R2
## 0.58102
## sla.R2_marginal.Marginal R2
## 0.16973
## gs.R2_conditional.Conditional R2
## 0.46869
## gs.R2_marginal.Marginal R2
## 0.39313
## amax.R2_conditional.Conditional R2
## 0.94110
## amax.R2_marginal.Marginal R2
## 0.49635## r2 models
models_q2 %>%
map(., AIC) %>%
unlist()## nodule_mass_ratio nodule_weight nodule_count wue_nlme
## 98.016 -65.217 62.791 302.027
## wue_log pnue_nlme pnue_log n_area_nlme
## 124.712 1048.375 85.397 88.499
## n_area_log sla gs amax
## -2.973 -925.487 -159.109 483.009Anova tables
map(models_q2, ~anova.lme(.x, type = "marginal"))## $nodule_mass_ratio
## numDF denDF F-value p-value
## (Intercept) 1 35 0.2300 0.6345
## treatment 3 35 0.8484 0.4769
## init_height 1 35 1.2452 0.2721
##
## $nodule_weight
## numDF denDF F-value p-value
## (Intercept) 1 41 3.298 0.0767
## treatment 3 41 4.623 0.0071
## init_height 1 41 1.376 0.2476
##
## $nodule_count
## numDF denDF F-value p-value
## (Intercept) 1 41 27.527 <.0001
## treatment 3 41 11.330 <.0001
## init_height 1 41 5.359 0.0257
##
## $wue_nlme
## numDF denDF F-value p-value
## (Intercept) 1 102 121.59 <.0001
## nfixer 1 6 0.34 0.5798
## treatment 3 102 66.94 <.0001
## init_height 1 102 0.63 0.4307
## nfixer:treatment 3 102 49.70 <.0001
##
## $wue_log
## numDF denDF F-value p-value
## (Intercept) 1 102 36.36 <.0001
## nfixer 1 6 4.09 0.0897
## treatment 3 102 26.88 <.0001
## init_height 1 102 0.40 0.5264
## nfixer:treatment 3 102 4.67 0.0042
##
## $pnue_nlme
## numDF denDF F-value p-value
## (Intercept) 1 102 89.71 <.0001
## nfixer 1 6 0.41 0.5455
## treatment 3 102 16.00 <.0001
## init_height 1 102 7.83 0.0061
## nfixer:treatment 3 102 6.70 0.0004
##
## $pnue_log
## numDF denDF F-value p-value
## (Intercept) 1 102 947.5 <.0001
## nfixer 1 6 0.5 0.4992
## treatment 3 102 3.6 0.0156
## init_height 1 102 2.2 0.1425
## nfixer:treatment 3 102 1.8 0.1541
##
## $n_area_nlme
## numDF denDF F-value p-value
## (Intercept) 1 102 15.472 0.0002
## nfixer 1 6 11.956 0.0135
## treatment 3 102 20.985 <.0001
## init_height 1 102 1.109 0.2947
## nfixer:treatment 3 102 2.026 0.1149
##
## $n_area_log
## numDF denDF F-value p-value
## (Intercept) 1 102 0.606 0.4381
## nfixer 1 6 15.964 0.0072
## treatment 3 102 7.372 0.0002
## init_height 1 102 0.575 0.4501
## nfixer:treatment 3 102 2.385 0.0735
##
## $sla
## numDF denDF F-value p-value
## (Intercept) 1 102 118.25 <.0001
## nfixer 1 6 0.43 0.5342
## treatment 3 102 0.25 0.8618
## init_height 1 102 6.32 0.0135
## nfixer:treatment 3 102 1.80 0.1521
##
## $gs
## numDF denDF F-value p-value
## (Intercept) 1 102 8.874 0.0036
## nfixer 1 6 6.401 0.0447
## treatment 3 102 18.662 <.0001
## init_height 1 102 1.015 0.3161
## nfixer:treatment 3 102 1.698 0.1722
##
## $amax
## numDF denDF F-value p-value
## (Intercept) 1 102 12.364 0.0007
## nfixer 1 6 7.969 0.0302
## treatment 3 102 0.257 0.8562
## init_height 1 102 5.555 0.0203
## nfixer:treatment 3 102 7.384 0.0002Post-Hoc: Tukey’s test
Maximal photosynthesis
as_tibble(emmeans(models_q2$amax,
pairwise ~ treatment*nfixer,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable()| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| no_additions nonfixer - plus_nutrients nonfixer | -0.1287 | 0.5156 | 102 | -0.2495 | 1.0000 |
| no_additions nonfixer - plus_water nonfixer | -0.4310 | 0.5208 | 102 | -0.8277 | 0.9911 |
| no_additions nonfixer - plus_water_nutrients nonfixer | -0.0716 | 0.4988 | 102 | -0.1436 | 1.0000 |
| no_additions nonfixer - no_additions fixer | -9.0888 | 3.2195 | 6 | -2.8230 | 0.2419 |
| no_additions nonfixer - plus_nutrients fixer | -7.0199 | 3.2267 | 6 | -2.1756 | 0.4646 |
| no_additions nonfixer - plus_water fixer | -11.4230 | 3.2202 | 6 | -3.5473 | 0.1120 |
| no_additions nonfixer - plus_water_nutrients fixer | -9.7834 | 3.2227 | 6 | -3.0357 | 0.1931 |
| plus_nutrients nonfixer - plus_water nonfixer | -0.3024 | 0.5379 | 102 | -0.5621 | 0.9992 |
| plus_nutrients nonfixer - plus_water_nutrients nonfixer | 0.0570 | 0.5175 | 102 | 0.1102 | 1.0000 |
| plus_nutrients nonfixer - no_additions fixer | -8.9601 | 3.2232 | 6 | -2.7799 | 0.2532 |
| plus_nutrients nonfixer - plus_nutrients fixer | -6.8913 | 3.2296 | 6 | -2.1338 | 0.4829 |
| plus_nutrients nonfixer - plus_water fixer | -11.2943 | 3.2243 | 6 | -3.5029 | 0.1174 |
| plus_nutrients nonfixer - plus_water_nutrients fixer | -9.6548 | 3.2280 | 6 | -2.9910 | 0.2025 |
| plus_water nonfixer - plus_water_nutrients nonfixer | 0.3594 | 0.5209 | 102 | 0.6900 | 0.9971 |
| plus_water nonfixer - no_additions fixer | -8.6577 | 3.2230 | 6 | -2.6863 | 0.2792 |
| plus_water nonfixer - plus_nutrients fixer | -6.5889 | 3.2302 | 6 | -2.0398 | 0.5257 |
| plus_water nonfixer - plus_water fixer | -10.9919 | 3.2235 | 6 | -3.4099 | 0.1296 |
| plus_water nonfixer - plus_water_nutrients fixer | -9.3524 | 3.2260 | 6 | -2.8991 | 0.2232 |
| plus_water_nutrients nonfixer - no_additions fixer | -9.0171 | 3.2193 | 6 | -2.8009 | 0.2476 |
| plus_water_nutrients nonfixer - plus_nutrients fixer | -6.9483 | 3.2267 | 6 | -2.1534 | 0.4743 |
| plus_water_nutrients nonfixer - plus_water fixer | -11.3513 | 3.2198 | 6 | -3.5255 | 0.1146 |
| plus_water_nutrients nonfixer - plus_water_nutrients fixer | -9.7118 | 3.2220 | 6 | -3.0142 | 0.1975 |
| no_additions fixer - plus_nutrients fixer | 2.0688 | 0.7107 | 102 | 2.9108 | 0.0810 |
| no_additions fixer - plus_water fixer | -2.3342 | 0.6752 | 102 | -3.4573 | 0.0175 |
| no_additions fixer - plus_water_nutrients fixer | -0.6947 | 0.6816 | 102 | -1.0191 | 0.9705 |
| plus_nutrients fixer - plus_water fixer | -4.4030 | 0.7131 | 102 | -6.1741 | 0.0000 |
| plus_nutrients fixer - plus_water_nutrients fixer | -2.7635 | 0.7236 | 102 | -3.8192 | 0.0055 |
| plus_water fixer - plus_water_nutrients fixer | 1.6395 | 0.6758 | 102 | 2.4259 | 0.2403 |
# Treatment effects
emmeans_table_tidy(models_q2$amax,
formula = "treatment|nfixer",
grouping_var = "nfixer")## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment | nfixer
## <environment: 0x5f61ade21730>Water Use Efficiency
as_tibble(emmeans(models_q2$wue_log,
pairwise ~ treatment*nfixer,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable() | contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| no_additions nonfixer - plus_nutrients nonfixer | 0.0227 | 0.1046 | 102 | 0.2167 | 1.0000 |
| no_additions nonfixer - plus_water nonfixer | 0.6581 | 0.1059 | 102 | 6.2117 | 0.0000 |
| no_additions nonfixer - plus_water_nutrients nonfixer | 0.6950 | 0.1015 | 102 | 6.8499 | 0.0000 |
| no_additions nonfixer - no_additions fixer | -0.4350 | 0.2152 | 6 | -2.0215 | 0.5342 |
| no_additions nonfixer - plus_nutrients fixer | -0.6114 | 0.2197 | 6 | -2.7829 | 0.2524 |
| no_additions nonfixer - plus_water fixer | -0.3098 | 0.2154 | 6 | -1.4384 | 0.8145 |
| no_additions nonfixer - plus_water_nutrients fixer | -0.2779 | 0.2164 | 6 | -1.2845 | 0.8769 |
| plus_nutrients nonfixer - plus_water nonfixer | 0.6354 | 0.1092 | 102 | 5.8190 | 0.0000 |
| plus_nutrients nonfixer - plus_water_nutrients nonfixer | 0.6723 | 0.1049 | 102 | 6.4101 | 0.0000 |
| plus_nutrients nonfixer - no_additions fixer | -0.4577 | 0.2171 | 6 | -2.1083 | 0.4943 |
| plus_nutrients nonfixer - plus_nutrients fixer | -0.6340 | 0.2213 | 6 | -2.8656 | 0.2313 |
| plus_nutrients nonfixer - plus_water fixer | -0.3325 | 0.2175 | 6 | -1.5291 | 0.7735 |
| plus_nutrients nonfixer - plus_water_nutrients fixer | -0.3006 | 0.2189 | 6 | -1.3736 | 0.8421 |
| plus_water nonfixer - plus_water_nutrients nonfixer | 0.0369 | 0.1059 | 102 | 0.3484 | 1.0000 |
| plus_water nonfixer - no_additions fixer | -1.0931 | 0.2173 | 6 | -5.0308 | 0.0254 |
| plus_water nonfixer - plus_nutrients fixer | -1.2694 | 0.2218 | 6 | -5.7235 | 0.0136 |
| plus_water nonfixer - plus_water fixer | -0.9679 | 0.2175 | 6 | -4.4509 | 0.0443 |
| plus_water nonfixer - plus_water_nutrients fixer | -0.9360 | 0.2183 | 6 | -4.2867 | 0.0522 |
| plus_water_nutrients nonfixer - no_additions fixer | -1.1300 | 0.2151 | 6 | -5.2535 | 0.0207 |
| plus_water_nutrients nonfixer - plus_nutrients fixer | -1.3063 | 0.2197 | 6 | -5.9465 | 0.0113 |
| plus_water_nutrients nonfixer - plus_water fixer | -1.0048 | 0.2153 | 6 | -4.6678 | 0.0358 |
| plus_water_nutrients nonfixer - plus_water_nutrients fixer | -0.9729 | 0.2161 | 6 | -4.5018 | 0.0421 |
| no_additions fixer - plus_nutrients fixer | -0.1764 | 0.1445 | 102 | -1.2206 | 0.9240 |
| no_additions fixer - plus_water fixer | 0.1252 | 0.1372 | 102 | 0.9120 | 0.9843 |
| no_additions fixer - plus_water_nutrients fixer | 0.1571 | 0.1380 | 102 | 1.1379 | 0.9468 |
| plus_nutrients fixer - plus_water fixer | 0.3015 | 0.1448 | 102 | 2.0828 | 0.4330 |
| plus_nutrients fixer - plus_water_nutrients fixer | 0.3334 | 0.1461 | 102 | 2.2824 | 0.3135 |
| plus_water fixer - plus_water_nutrients fixer | 0.0319 | 0.1372 | 102 | 0.2325 | 1.0000 |
# Treatment effects
emmeans_table_tidy(models_q2$wue_log,
formula = "treatment|nfixer",
grouping_var = "nfixer")## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment | nfixer
## <environment: 0x5f61aba70700>as_tibble(emmeans(models_q2$wue_log,
pairwise ~ nfixer,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable()| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| nonfixer - fixer | -0.7525 | 0.189 | 6 | -3.981 | 0.0073 |
# Treatment effects
as.data.frame(emmeans::emmeans(models_q2$wue_log,
specs = pairwise ~nfixer,
type = "response",
adjust = "tukey")$emmeans) %>%
janitor::clean_names() %>%
dplyr::select(response, everything(),
# Remove variables
-c(df, lower_cl, upper_cl, se)) %>%
# Rename response to emmean, this is done when models is log
dplyr::rename_all(funs(stringr::str_replace_all(., "response", "emmean"))) %>%
# Calculate % difference between control and variable, this assume that
# that first name is the control
dplyr::mutate(difference = ((emmean - first(emmean))),
perc_difference =((emmean - first(emmean) )/first(emmean))*100) %>%
dplyr::mutate_if(is.numeric, round, 3)## emmean nfixer difference perc_difference
## 1 1.888 nonfixer 0.000 0.0
## 2 4.008 fixer 2.119 112.2Stomatal Conductance
as_tibble(emmeans(models_q2$gs,
pairwise ~ treatment,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable() | contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| no_additions - plus_nutrients | 0.0282 | 0.0247 | 102 | 1.1424 | 0.6643 |
| no_additions - plus_water | -0.1245 | 0.0239 | 102 | -5.1985 | 0.0000 |
| no_additions - plus_water_nutrients | -0.1041 | 0.0236 | 102 | -4.4076 | 0.0002 |
| plus_nutrients - plus_water | -0.1527 | 0.0251 | 102 | -6.0899 | 0.0000 |
| plus_nutrients - plus_water_nutrients | -0.1323 | 0.0248 | 102 | -5.3291 | 0.0000 |
| plus_water - plus_water_nutrients | 0.0204 | 0.0239 | 102 | 0.8507 | 0.8300 |
# Treatment effects
emmeans_table_tidy(models_q2$gs,
formula = "treatment",
)## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x5f61a605f958>SLA
SLA did not responded to any of the treatments nor varied between nfixers and non-fixing species so no post hoc test was performed.
PNUE
as_tibble(emmeans(models_q2$pnue_log,
pairwise ~ treatment,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable()| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| no_additions - plus_nutrients | 0.2100 | 0.0745 | 102 | 2.8193 | 0.0291 |
| no_additions - plus_water | -0.1149 | 0.0722 | 102 | -1.5899 | 0.3889 |
| no_additions - plus_water_nutrients | 0.0079 | 0.0715 | 102 | 0.1111 | 0.9995 |
| plus_nutrients - plus_water | -0.3248 | 0.0758 | 102 | -4.2844 | 0.0002 |
| plus_nutrients - plus_water_nutrients | -0.2020 | 0.0755 | 102 | -2.6757 | 0.0425 |
| plus_water - plus_water_nutrients | 0.1228 | 0.0722 | 102 | 1.6998 | 0.3291 |
# Treatment effects
emmeans_table_tidy(models_q2$pnue_log,
formula = "treatment",
)## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x5f61ad643908>Nitrogen concentration per unit of area
as_tibble(emmeans(models_q2$n_area_log,
pairwise ~ treatment,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable()| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| no_additions - plus_nutrients | -0.1575 | 0.0480 | 102 | -3.2823 | 0.0076 |
| no_additions - plus_water | -0.0236 | 0.0465 | 102 | -0.5074 | 0.9572 |
| no_additions - plus_water_nutrients | -0.0502 | 0.0461 | 102 | -1.0895 | 0.6968 |
| plus_nutrients - plus_water | 0.1339 | 0.0489 | 102 | 2.7383 | 0.0361 |
| plus_nutrients - plus_water_nutrients | 0.1073 | 0.0489 | 102 | 2.1925 | 0.1322 |
| plus_water - plus_water_nutrients | -0.0266 | 0.0465 | 102 | -0.5722 | 0.9402 |
# Treatment effects
emmeans_table_tidy(models_q2$n_area_log,
formula = "treatment",
)## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x5f61b0b25188>as_tibble(emmeans(models_q2$n_area_log,
pairwise ~ nfixer,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable()| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| nonfixer - fixer | -0.7214 | 0.1866 | 6 | -3.865 | 0.0083 |
# Treatment effects
as.data.frame(emmeans::emmeans(models_q2$n_area_log,
specs = pairwise ~nfixer,
type = "response",
adjust = "tukey")$emmeans) %>%
janitor::clean_names() %>%
dplyr::select(response, everything(),
# Remove variables
-c(df, lower_cl, upper_cl, se)) %>%
# Rename response to emmean, this is done when models is log
dplyr::rename_all(funs(stringr::str_replace_all(., "response", "emmean"))) %>%
# Calculate % difference between control and variable, this assume that
# that first name is the control
dplyr::mutate(difference = ((emmean - first(emmean))),
perc_difference =((emmean - first(emmean) )/first(emmean))*100) %>%
dplyr::mutate_if(is.numeric, round, 3)## emmean nfixer difference perc_difference
## 1 1.037 nonfixer 0.000 0.0
## 2 2.133 fixer 1.096 105.7Nodule weight
as_tibble(emmeans(models_q2$nodule_weight,
pairwise ~ treatment,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable()| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| no_additions - plus_nutrients | -0.0631 | 0.0231 | 41 | -2.7291 | 0.0443 |
| no_additions - plus_water | -0.0227 | 0.0129 | 41 | -1.7523 | 0.3107 |
| no_additions - plus_water_nutrients | -0.0605 | 0.0198 | 41 | -3.0512 | 0.0200 |
| plus_nutrients - plus_water | 0.0404 | 0.0239 | 41 | 1.6944 | 0.3397 |
| plus_nutrients - plus_water_nutrients | 0.0026 | 0.0265 | 41 | 0.0978 | 0.9997 |
| plus_water - plus_water_nutrients | -0.0378 | 0.0205 | 41 | -1.8480 | 0.2663 |
# Treatment effects
emmeans_table_tidy(models_q2$nodule_weight,
formula = "treatment")## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x5f61ad9f9038>Nodule Count
as_tibble(emmeans(models_q2$nodule_count,
pairwise ~ treatment,
adjust = "tukey"
)$contrast) %>%
mutate(across(2:6, round, 6)) %>%
kable()| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| no_additions - plus_nutrients | 0.1686 | 0.0642 | 41 | 2.6250 | 0.0565 |
| no_additions - plus_water | 0.1165 | 0.0689 | 41 | 1.6909 | 0.3415 |
| no_additions - plus_water_nutrients | -0.1858 | 0.0630 | 41 | -2.9504 | 0.0258 |
| plus_nutrients - plus_water | -0.0522 | 0.0648 | 41 | -0.8059 | 0.8513 |
| plus_nutrients - plus_water_nutrients | -0.3545 | 0.0640 | 41 | -5.5420 | 0.0000 |
| plus_water - plus_water_nutrients | -0.3023 | 0.0683 | 41 | -4.4277 | 0.0004 |
# Treatment effects
emmeans_table_tidy(models_q2$nodule_count,
formula = "treatment")## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x5f61b0a5f530>Nodule mass ratio
Nodule mass ratio did not responded to any of the treatments applied so no post hoc test was performed.
Boxplots
Boxplots traits
# Step done for getting predictions from models for Q2
data_for_predictions <-
data_for_models %>%
rownames_to_column("id") %>%
# Remove unused variables
dplyr::select(id, spcode, treatment, nfixer, init_height)string <- c("models_q2")
data_pred_traits <-
# Get models prediction
gather_predictions(data_for_predictions ,
# Return predictions
models_q2$amax,
models_q2$sla,
models_q2$gs,
models_q2$wue_log,
models_q2$pnue_log,
models_q2$n_area_log
) %>%
pivot_wider(names_from = model, values_from = pred) %>%
rename_all(funs(
# rename columns
stringr::str_to_lower(.) %>%
stringr::str_replace(., c(string),"pred_") %>%
# Remove dollar sing
gsub("\\$", "", .)
)) %>%
# Back transform log variables
mutate(pred_wue = exp(pred_wue_log),
pred_n_area = exp(pred_n_area_log),
pred_pnue = exp(pred_pnue_log),
) %>%
# Remove log predictions and init height
dplyr::select(-c(init_height, pred_wue_log, pred_n_area_log))# Generate plot combinations
vars_q2_interaction <-
crossing(
# Get all numeric variables to plot (all y)
as_tibble(t(combn(dplyr::select(data_pred_traits, where(is.numeric)) %>% names, 1))),
# Select factor variables to plot
x_axis_var = dplyr::select(data_pred_traits, nfixer) %>% names,
group_var = dplyr::select(data_pred_traits, treatment) %>% names) %>%
filter(V1 %in% c('pred_amax', 'pred_wue'))vars_q2_interaction %>%
# Gererate plots
pmap( ~ boxplot_plot_pmap(data = data_pred_traits,
y = !!sym(..1), x = !!sym(..2),
fill = !!sym(..3)))## [1] 24.08
## [1] 7.034## [[1]]
##
## [[2]]
vars_q2_treatment <-
crossing(
# Get all numeric variables to plot (all y)
as_tibble(t(combn(dplyr::select(data_pred_traits, where(is.numeric)) %>% names, 1))),
# Select factor variables to plot
x_axis_var = dplyr::select(data_pred_traits, treatment) %>% names,
group_var = dplyr::select(data_pred_traits, treatment) %>% names) %>%
filter(V1 %in% c('pred_gs', 'pred_n_area', 'pred_pnue' ))vars_q2_treatment %>%
# Gererate plots
pmap( ~ boxplot_plot_pmap(data = data_pred_traits,
y = !!sym(..1), x = !!sym(..2),
fill = !!sym(..3)))## [1] 0.359
## [1] 3.399
## [1] 117## [[1]]
##
## [[2]]
##
## [[3]]
Boxplots Nodules
data_for_nodule_predictions <-
data_nodules_cleaned %>%
rownames_to_column("id") %>%
# Remove unused variables
dplyr::select(id, spcode, treatment, init_height)string <- c("models_q2")
data_pred_nodules <-
# Get models prediction
gather_predictions(data_for_nodule_predictions,
# Return predictions
models_q2$nodule_count,
models_q2$nodule_mass_ratio,
models_q2$nodule_weight) %>%
pivot_wider(names_from = model, values_from = pred) %>%
rename_all(funs(
# rename columns
stringr::str_to_lower(.) %>%
stringr::str_replace(., c(string),"pred_") %>%
# Remove dollar sing
gsub("\\$", "", .)
)) %>%
# Back transform log variables
mutate(pred_nodule_count = exp(pred_nodule_count),
pred_nodule_mass_ratio = exp(pred_nodule_mass_ratio)) %>%
# Remove log predictions and init height
dplyr::select(-c(init_height))# Generate plot combinations
vars_q2_nodules <-
crossing(
# Get all numeric variables to plot (all y)
as_tibble(t(combn(dplyr::select(data_pred_nodules, where(is.numeric)) %>% names, 1))),
# Select factor variables to plot
x_axis_var = dplyr::select(data_pred_nodules, treatment) %>% names,
group_var = dplyr::select(data_pred_nodules, treatment) %>% names) vars_q2_nodules %>%
# Gererate plots
pmap( ~ boxplot_plot_pmap(data = data_pred_nodules,
y = !!sym(..1), x = !!sym(..2),
fill = !!sym(..3)))## [1] 167.2
## [1] 2.618
## [1] 0.199## [[1]]
##
## [[2]]
##
## [[3]]